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1.  INTRODUCTION 


A  decade  ago  some  facts  and  fallacies  in  digital  simulation  were  discussed  in  an 
interesting  article  [  1  ].  Some  of  the  facts  appear  to  be  fallacies  as  that  author  had 
warned  might  be  the  case.  The  fallacies  arose,  in  part,  because  of  premature  use  of 
the  general  recurrence;  that  is,  improper  incorporation  of  the  initial  conditions.  It. 
is,  of  course,  unfair  to  “sharp  shoot "  after  a  decade,  but  some  of  the  “facts"  (fallacies) 
are  still  present  in  today’s  z-transform  literature. 

These  and  other  questions  were  addressed  in  an  earlier  report  [2]  documenting 
the  results  of  an  initial  investigation. 

Several  suggestions  were  made  as  to  avenues  for  further  possible  development. 
They  were  as  follows: 

a)  Place  z-transforms  on  a  firm  foundation  using  distribution  theory, 

b)  Determine  when  modified  z-transforms  and/or  tunable  convolution  arc 
advantageous  and  in  what  combination.  This  would  include  the  use  of  higher  order 
holds. 

c)  Analyze  the  effects  of  tuning  for  other  inputs  and  other  transfer  functions. 

Some  progress  has  been  made,  and,  though  in  no  way  complete  at  this  time,  it 
was  felt  that  the  results  to  date  were  worth  reporting  if  only  to  stimulate  further 
interest. 

There  is  the  story  of  a  physicist  and  a  mathematician  in  an  airplane  which  flew 
over  a  flock  of  sheep,  all  white  save  one,  who  was  black.  The  physicist  proceeded  to 
theorize  about  the  number  of  black  sheep  in  the  universe.  The  mathematician  knew 
there  was  one  flock  with  a  sheep  who  was  black  ON  TOP.  An  engineer  would 
probably  choose  to  ignore  the  color  of  the  sheep  since  it  would  have  no  effect  on  the 
mutton.  This  author  is  a  physicist. 

2.  SOME  PRELIMINARIES 

In  this  section  certain  details  as  to  the  notation  used  and  some  conjectures 
concerning  the  use  of  distributions  in  z-transforms  will  be  presented. 
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Thai  ihc  l.aplace  transform  of  the  Dirac  delta*  is  a  constant  and  the  Laplace 
transform  of  the  n-th  derivative  of  the  Dirac  delta  is  s"  (see  Appendixes  A  and  B) 
would  certainly  seem  to  indicate  an  intimate  connection  between  Laplace  transforms 
and  distribution  theory. 

For  reasons  which  will  become  apparent,  the  Laplace  transform  will  be  defined  as 
13]. 


<v 

L(f(t)>  /  f (t >  u(t)  e~8t  dt 

«U) 


(1) 


where  u(t)  is  the  Heaviside  unit  step  (whose  derivative  is  the  Dirac  delta). 
Of  course  Equation  (1)  readily  reduces  to 


oo 

f(s)  -  /  f  (t )  e~stdt 

o 


(2) 


Note  that  the  introduction  of  the  unit  step,  Equation  (1),  remov  i  the  need  for 
defining  the  behavior  of  “fU>”  prior  to  “t”  equal  zero.  The  lower  lim.  of  integration 
in  Equation  (2)  is  contained  in  the  unit  step  in  Equation  (1).  For  example,  consider 
convolution: 


OO 


g(t)  *  f (t)=  /  g(t) 

.oo 

u(t)  f(t-t)  u(i-t)  dt. 

(3) 

oo 

-  /  g(t) 

f(t-t)  u(t)  u(T-t)  dt. 

(4) 

.OO 


The  term,  “u(t)  u(r— t),”  will  be  recognized  as  the  unit  pulse  from  0  to  r. 
Equation  (4)  becomes 

i 

g(t)  *  f(t)  -  /  g(t)  f(T-t)  dt. 
o 


•The  Dirac  delta  distribution  is  not  a  function.  Dirac  to  distinguish  it  from  the  Kroneker  delta. 


(5) 


6 


It  should  be  noted  that  the  Dirac  delta  is  the  unity  element  in  convolution, 

on 

/  g(t)  6(i-t)  dt  «  g(0.  (6) 

»  a> 

For  a  discussion  of  the  convolution  of  distributions,  see  Kecs  and  Teodorescu  [4], 
The  Dirac  delta  may  also  be  thought  of  as  a  sampler 

Oil 

/  g(t)  6(nT-t )  dt  =  g(nT)  (7) 

..On 

If  one  has  the  sequence  [5] 

on 

f*(t)  =  [  f(t)  6(nT-t)  (8) 

n-0 

its  Laplace  transform  is 

OO 

L|f*(t)l  *  l  f(nT)  e  nsT  •  (9) 

n«0 

Taking 

z  =  e-st  (10) 


one  has 

L  [  f *Ct )  ]  *  l  f(nT)  zn  .  (11) 

n»0 
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The  z-transform  may  be  thought  of  as  a  discrete  Laplace  transform,  that  is.  a 
transformation  to  the  sequency  domain.  The  notation 


z(f(8)l  L[  f*(Oi  (12) 

will  be  used  in  the  following.  To  be  consistent,  let 

Z  [ f  ( s )  1  =  [  f(nT)  ti  ( nT )  2"  (13) 

n=~^ 

and 

oi' 

zlf(t))  =  l  f  (t)  ,S(nT-t)  •  (14) 

n=— x’ 

Since 

g(s)  t (s)  =  Llg(t)  *  f (t)J  (15) 

it  follows  that 

Z[g(s)  f(z)l  = 

00  CO  00 

£  z"  u(nT)  /  (g(nT-t)  u(nT-t)  £  f(kT)  u(t)  6(kT-t)J  dt 
n=-'x>  -<»  k5*-* 

(16) 

From  the  properties  of  the  Dirac  delta,  [6,7]  and  Appendix  B, 
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(17) 


ll 

n  ■ 


*  I  *  g (nT-kT)  f(kT) 

n*o  k-o 


l  g(ni-kT)  /”  k  1  (kT)  7.k 

n»0  k”0 


(18) 


which  is  the  Cauchy  product  of  power  series 


)  ;  t 

j  \  g(i»T)  zn,|}.  t  (kT)  zkj 
(  n*0  j  \k  =  0 


(19) 


su)  t  ( z ) , 


(20) 


the  Ragazzini-Zadeh  identity. 

Zlg(s)  i(s)]«  y ( 2 ) t ( z ) 

but  unfortunately 

*•  [  g t S )  f  ( s )  ]  g  (  z  )  f(z) 

that  is 

gU)  #  y(z), 

and  therein  lies  the  difficulty  in  applying  z-transforms  to  continuous  systems,  that 
is,  the  digital  simulation  of  continuous  systems. 

3.  MEAN  VALUE  CONVOLUTION 


A  very  useful  relationship  would  be  the  solution  of 

OO  CD 

Z(g(s)  f(s>)  -  l  zn  u(nT)  /  g(t)  u(t)  f(nT-t)  u(nT-t)  dt  (21) 

—CD 
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where  both  “g(t)”  and  are  functions.  There  are  no  Dirac  deltas  to  “remove” 
the  integral  as  was  the  case  with  the  Ragazzini-Zadeh  identity. 

The  mean  value  theorem  of  the  integral  calculus  [5]  guarantees  there  is  some  “6k” 
such  that 

(k+l)T 

/  g(t)  f(nT-t)  dt  =  T  g(kT  +  6  T)  f(nT-kT-6kT)  (22) 

kT 

and  Equation  (21)  may  be  rewritten 

oo  n-1 

l  zn  l  T  g(kT+6  T)  f (nT-kT-6,  T)  (23) 

n=0  k=0 


Since  “6k”  is  most  likely  different  for  each  “k”  it  would  be  difficult  to  proceed. 
Assuming 


6k  =  6  (24) 


that  is,  “6k”  is  the  same  for  all  “k”  one  has 


oo  n-1 

l  zn  l  T  g(kT-<5T)  f(nT-kT-6T)  (25) 

n=0  k=0 


and  as  before 


oo  Il“l 

T l  l  g(kT+6T)  zk  f(nT-kT-6T)  zn_k  *  (26) 

n»0  k=0 

One  problem,  the  sum  is  to  "n  -  1”  not  ”n”and  to  obtain  the  form  of  the  Cauchy 
product  a  term  must  be  added  and  subtracted. 
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I  1 

V  K(kT+iST)  zR 

I}'  1  (nT-iST)  z‘* 

-f (-6T) 

k*0 

[n-0 

and.  finally,  from  the  definition  of  the  modified  z-transform 

Z|n(s)  t(s)|  T  g(z.vS)|f(z,-6)  -  t  (-6T)  1  (28) 

which  will  be  referred  to  as  Mean  Value  Convolution,  MVC  for  short. 


Since  the  assumption.  Equation  (24),  is  suspect,  some  checks  arc  warranted.  > 
Consider  the  case  where 


(29) 


and 


From  Equation  (28)  and  Appendix  A 


(30) 


(31) 


Checking  against  Appendix  A  for  I/s2  Equation  (31)  is  found  to  be  exact.  Forg(s)the 
same  and 
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(32) 


f(s) 


t 


2 

s 


Equation  (28)  and  Appendix  A  yield 


i'4  ter 


nr 


) 


which  simplifies  to 


.I"2  (( ten)  -  nz  I 
(  1  -  2)  3 

For  a  small  time  step,  one  would  suspect  that 


(33) 


(34) 


(35) 


would  be  reasonable  and  since 


n  -  -  5 


one  would  have 

T2z  (l+z) 
2(I-z)3 

which  for  I/s'  is  exact. 
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(36) 


(37) 


Interchanging  the  roies  of  g(s)  and  fls)  in  Equation  <28),  one  has 


and  simplifying 

T2z  (n  +  (l  -  T))  2) 
<1- z)3 


(38) 


(39) 


In  this  case 


n  -  6  -  y  (40) 

and  Equation  (39)  becomes 

72*XL+Z)  (41) 

2(1  -  z)3 

which  for  I/s1  is  again  exact. 

That  the  above  checks  are  exact  is  not  surprising  since  the  first  would  correspond 
to  integration  of  a  constant,  and  the  second  to  integration  of  a  ramp  or  the  double 
integration  of  a  constant.  For  these  cases 


6^  =  6  *  (42) 

is  what  one  would  expect. 

So  far,  so  good,  but  l/s4  presents  difficulties  whose  discussion  will  be  deferred  until 
a  solution  to  the  difficulty  is  available  in  a  later  section. 
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4.  R-K  CONVOLUTION 


In  the  previous  section  the  mean  value  theorem  was  invoked  to  “remove"  the 
integral.  In  this  section  the  approach  will  be  to  “numerically"  integrate.  Since  the 
approach  of  this  report  is  that  of  z-transforms  a  brief  summary  of  integrators  with 
holds  is  in  order.  Consult  Jury  [6),  Smith  [5]  or  Rosko  |7|  for  details. 


Z[s(s)  f  (s)]  -  y(z)  f  (z) 


(43) 


may  be  approximated  by 


Z Qi ( s )  t  (s)] 


g(z)  -  S 

7. 


t&) 


(44) 


that  is 


y(*!>  -  g(z) 


t 


lH-1  f 


f  (z)  Z 


(M 


(45) 


where  “n”  is  the  order  of  the  hold. 


For  an  integrator,  “f(s)  =  1/s",  and  a  zero  order  hold,  “n  =  0",  one  has 


«(z) 


z(?)  &)8U> 
z(»)  (^) 


1  z  JSJ& 

1  -  z 


(46) 
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Eular  integration  or  the  Left  Riemann  Sum. 
For  a  first  order  hold,  “n  =  1  ”, 


TO  +z) 
2(1  -z) 


g(z)  . 


(47) 


trapezoidal  integration. 

For  a  second  order  hold,  “n  =  2”, 

/t3(1  Mz  Iz2 
V  6(1  -  z)4 

/  T2z(1  4  z)\ 

\2(l-z)3  / 

Simpson's  rule. 

Smith’s  “tunable”  integrators  [5]  may  be  obtained  by  using  the  modified  z~ 
transform.  The  “tunable"  integrators  for  zero  and  first  order  holds  will  be  found  in 
Table  /. 

One  may  proceed  to  higher  order  holds  ( Table  1),  but  these  integrators  do  not 
appear  practical. 


, 


l 

g'( 


T  1  +  4z  + 


3  (1  +  z)  (l~z) 


g(z) 


(48) 


A  digression: 

Of  course,  holds  may  be  used  with  plants  other  than  integrators;  for  example,  a 
single  pole  filter, 


f(s> 


1 

8+a 
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(49) 


A  single  pole  filter  and  a  zero  order  hold,  "n  =  0”,  yield 

z(s  (s  +  a))  _  ( 1  -•£»  a^)  z  k(z) 


g(z) 


a(l-z  e"aT) 
If  one  chooses  to  implement  Figure  1, 


(50) 


g(s)  - 


Figure  1.  An  Integrator  with  feedback. 

an  integrator  with  feedback  for  a  zero  order  hold,  one  would  have 

z/g(s)\  T  z  g(z) 

\s  )  l-z(l-aT) 

which  would  amount  to  using  a  (1/0)  Pade'  approximation  for  e"*T, 


(51) 


(52) 


in  Equation  (50). 

Dividing  Equation  (51)  by  (50)  one  has  the  ratio* 


/  aT  \l 

'i  -  2e-T  \ 

\l -e_aT/ 

^l-z(l-aT)/ 

'Subtract  on*  and  multiply  by  on*  hundred  for  percent  error. 
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Applying  the  final  value  theorem,  /  —I  (n  the  ratio  is  one.  as  desired,  but 

applying  the  initial  value  theorem.  /  —0  (n— 0).  one  has  the  ratio 


Since  “a”  is  given,  **T”  must  be  chosen  so  that  Equation  (52)  is  as  near  one  as 
desired,  or  use  Equation  (50)  instead. 


For  a  single  pole  filter  and  a  first  order  hold  one  has 

*(  v  1  ) 

\s_(s  +  a)/ 

w 


«u> 


-aT  -aT  -aT 

,  lAuT+e  .  .  -JA  +  ( 1-e  -ate.  .  >zJ&U> 


>  — aT 

a‘T(l-n>  ) 


(53) 


It  is  common  practice  in  digital  simulation  to  follow  the  analog  computer 
practice  of  reducing  problems  to  interconnected  single  integrator  \  This  is  not 
necessary  nor  may  it  be  desirable,  but  enough  of  this  digression. 

In  the  following, several  integrators  (Table  1)  will  be  used  to  find  i.oproximate 
solutions  to  Equation  (21).  There  is  some  repetition  of  earlier  material  |21  which  is 
included  for  completeness. 

Using  Eular  integration.  Equation  (46),  in  Equation  (21)  one  has* 

n-l  (k  +  1)^  ti-l 

y  Zn  V  /  g(T)  t(nT-r)  dt  =  y  y  T  g(kT>  f(nT-kT)  .  (54) 

n»0  k“0  kT  n=0  k=0 


As  before,  adding  and  substracting  .  .  . 


oo 


-  I  T 

n*0 


R(kT) 


t.k-0 


f  (nT-kT)-g(nT)  fQ  . 


(55) 


*How  recurrence*  nre  developed  i»  treated  in  the  sample  problem  section. 
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The  Cauchy  product .  . 


OO 

OO 

T 

l  g(kT)  zk 

l  f (nT)  Zn 

k=0 

,n=  0 

-TfJ  g(nT)  zn  , 


n=0 


and  from  the  definition  of  the  z-transform, 


(56) 


Z[g(s)f(s)]=  T  g(z)  [f(z)-fo] 


(57) 


Eular  convolution. 

If  the  modified  z-transform  is  used  [2],  for  a  zero  order  hold  one  has 

°°  n  n_1 

l  z  Tl  (n  g(nT+T)  f(nT-kT-T)  +  (1-n)  g(nT)  f(nT-kT)] 
n=0  k=0 


(58) 


and  then,  as  above. 


oo  n 


-  tL 


n=0  k=0  |  °  n*0  °  n=0 


T l  zn  l  g(kT)  f(nT-kT)  -  T|n  gn  I  z"  f  (nT)  +  (1-n) fft  l  z"  g(nT)J 

(59) 


and  finally, 


zfe(s)  f(sj]  =  T  g(z)  f(z)-T[ng0  f(z)  +  (l-n)  g(z)  f  ]  ,  (60) 
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Eular  tunable  convolution. 


Proceeding  on  with  the  trapezoidal  integrator, 


n-l 


n  T 

1  2 

n-0  k-0 


l  zn  \  l  [g(kT+T)  f (nT-kT-t)  +  g(kT)  f(nT-kT) 


(61) 


and  after  manipulating  indicies 


l  znT  l  g(kT)  f(nT-kT)  -  y(g  f(nT)  +  g(nT)  f  ] 
n-0  k-0  0 


(62) 


and  rearranging 


00  n  .  , 

TV  V  g(kT)  z  f(nT-kT)  zn" 
n-0  k-0 


1 

2 


[f  l 

[  °  n-0 


g (nT)  +  g  I  zn  f (nT) 


n-0 


(63) 


and  finally . 


Z[g(s)f(s)]  =  T  f(z)  g(z)  -  f o  g(z)  +  gQ  f(z)] 


(64) 


trapezoidal  convolution  [7,8,9].  Equation  (64)  may  also  be  written 


z[g(s)  f(s)]=  y  [(g(z)  -  g  )  f(z)  +g(z)  (f(z)  -  fc)] 


(65) 
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For  the  modified  z-transform  and  first  order  hold  ( Table  1 )  one  has 

'•»•>  n-1  r 

l  Zn  |  )'  n2g(kT+2T)  f(nl-kT-2T) 

n=()  k  =  ()L 

+  (l  +  2n-2n2)  g (kT+T)  f(nT-kT-T)  +  (1-2  n  +  n2)  g(kt)  f(nT-kT)j 

(66) 

Changing  summing  indices 

/  ,0+1  ‘  n 

iV  znJrfjl  g(kT)  f(nT-kT)  +  (l+2n-2if)£  g(kT)  f(nT-kT) 

•  n-0  I  k=2  k-1 


+  (l-2n  +  n2)I  g(kT)  f(nT-kT)?  (67) 
k-0  / 

and  adding  and  subtracting  the  necessary  terms, 

OL’  l  n 

t!'  l  g(kT)  f(nT-kT)  -  |(1  +  2n  +  n2)  g0  f(nT) 

n=0  (  k»0 

-  -|(1  -  2n  +  n2)  fQ  g(nT) 

+  ~  n2  [g (nT+T)  f(-T)  -  g(T)  f (nT-T)l 

(68) 
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From  the  definition  of  the  z-transform  (and  the  Cauchy  product),  finally, 


Z[g(s)t(s)]  t  *u)  t  (.’)  -'■(!+  2n  -  .f)g  f  (z) 

-  id  -  2*.  +  -i2)f  f(z) 

2  o 

T  •  [«(/,!>  f  ( -T )  -  g(T)  f  (z,-l )]  .  (69) 

’T  T 

Trapezoidal  tunable  convolution. 

Simpson’s  rule  (48), 

T  1  +  l*  z  -»-  z*~  , 

3  (1  +z)  (1  -  Z) 


presents  indexing  problems  since  the  recurrence  would  be 


X 

n 


X  „  +  ?[ X  + 
n-2  3  n 


4Xn-l  +  Xn-21 


(70) 


Also,  Simpson’s  rule  is  known  to  be  sensitive  to  noise  [6,10].  Another  form,  used  in 
third  order  Runge-Kutta  integration,  is 


T  1  +  4z1/2  +  z 
^  (1  -  z) 


(71) 


whose  recurrence  is 


n-1 


1  fX 
6  1  n 


+  4X 


-1/2 


Xn-1] 


(72) 
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Equation  (71)  may  be  factored  into 


1  llLlf  \  +  2  /t  z1/2\ 

3  [2(1  -  zjj  3  -  z  y 


that  is. 


>'n  =  2'n-l  +  1  [7(>n+>'n-l)j+  f  T  t  n-1/2 


Note  that  Equations  (73)  and  (74)  are  a  linear  combination  of  trapezoidal 
integration,  and  Mean  Value  integration  with  <5  =  I  2.  Therefore  R-K  convolution  is  a 
linear  combination  of  trapezoidal  convolution  and  Mean  Value  Convolution, 


Z  jg(s)  its)  ]  =  |  j(g(z)  - 


£u)  f(z) 


+  t.g(a,l/2)  If  (z,-l/2)-  f  (-T/2 )] 


+  g(z)  (f (z)  -  f 


5.  ANALYSIS  OF  ERRORS 

Earlier  the  ratio  between  the  zero  order  hold  integrator  with  feedback  and  the 
zero  order  hold  single  pole  filter  was  taken  to  compare  these  different  simulation 
approaches.  Equation  (51).  The  input  was  not  specified. 

Now  the  ratio  of  the  convolution  approximation  to  the  exact  solution  will  be 
taken  for  given  inputs.  The  accuracy  is  highly  dependent  upon  the  input. 


Checks  for  the  various  convolution  approximations,  Table  2,  similar  to  those 
done  for  MVC,  Equations  (29)  through  (42),  will  be  found  in  Table  3.  The  ratios  to 
the  exact  z-transform  will  be  found  in  Table  4. 
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TABLE  2.  CONVOLUTION  APPROXIMATIONS  FOR  Z[g  (»)  t  (*)] 


T  fe(z,6)  [f (z,-6)  -  f (-6T)  ] 

T  g(z)  [f(z)  -  foJ 

|  [(g(z)  -2ngfl)  f  (z)  +  g(z)  (f(z)  -  2  (1-n)  fQ)] 

^  [(g(z)  -  gQ)  f(z)  +  g(z)  (f(z)  -  fQ)J 

^|f(g(z)  -  (i+2n-n2)  goJ  f(z)  +  g(z)£f(z)  -  (i-2n+n2) 

+  n2 Jg(z,l)  f(-T)-  g(T)  f(z,-l)]| 

I  j(g(z)-g0)  f(z)  +  4g(z,l/2)  [f(z,-l/2)-  f (-T/2^  +  ’(*>  (f(z)-f0) 


TABLE  3.  CHECKS  OF  THE  CONVOLUTION  APPROXIMATIONS  FOR  f(s)  =  '  AND 
g(t)  A  CONSTANT,  RAMP  OR  PARABOLA 


TABLE  4.  RATIOS  OF  APPROXIMATE  SOLUTION  TO  THE  EXACT  FOR  A 
CONSTANT,  RAMP  AND  PARABOLIC  INPUT  TO  A  SINGLE 
INTEGRATOR 


I  Ik-  ratios  tor  I  ulur  convolution  tevpme  some  interpretation  l  or  i  -  0  the  ratio 
would  he  zero  lor  the  tamp  ami  parabolic  inputs.  II  the  in  the  numerator  is 
interpreted  tin  ti  shift,  the  ratios  would  he  two  and  three,  respectively.  The 
approximation  for  the  ramp  was  determined  with 


r  i  r  l  I  .’ 

\-f  i  1  * 

|u-*)-Jli  '  J  "ti -<> * 


lnterehanKinK  K(s)  and  (\s)  in  Kquation  (f>7)  one  has 


Lu-*r  J  o  -*) 


and  the  ratio  would  be 


(76) 


(77) 


(»  +  e) 


(78) 


But  recalling  that 


T 

(1-2) 


(79) 


ia  Rectangular  integration,  one  may  conclude  the  difference  between  Eular 
integration  and  rectangular  integration,  the  Right  Riemann  Sum,  is  a  time  step. 
This  difficulty  arises  because  of  all  the  approximations  in  Table  2,  only  Eular 
Convolution  lacks  symmetry,  and  convolution  should  be  symmetric. 


One  could  proceed  to  higher  powers  of  time  for  inputs,  but  b  ing  only  able  to 
determine  the  ratio  at  the  initial  time  is  not  very  enlightening.  A  time  step  by  time 
step  approach  is  possible  |2|  but  does  not  seem  practical. 

Fortunately  some  transcendental  functions  work  very  well  as  inputs  for  ratio 
analysis.  For  example,  consider  an  exponential  input  to  an  inte,  rator.  That  is, 


f(s)  *  - 
8 

(80) 

...  -at 

g(t)  =  a  e 

(81) 

and,  it  follows  that 


g(s) 


a 

s  +  a  * 


(82) 


taking  the  z-transform  of  Equations  (80)  and  (82), 


f  (z) 


1 

1-2  ’ 


(83) 
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g(z) 


a 


1  -  z  t' 


-aT 


and  substituting  into  Equation  (57)  yields 


aT  z 

(1  -  z)  (1  -  z  e  aT) 


for  Eular  Convolution. 
Taking  the  z-transform  of 


7 /  a  \ _ ( 1  -  t‘  ~a  1  )z 

Vs(s+a)/  (l-z)(l-ze-aT) 

yields  the  exact  solution. 

Taking  the  ratio  of  Equation  (85)  and  (86), 


aT 

i  ~at 

1-e 

Note  that  no  “z’s”  remain. 

When  this  occurs,  the  assumption, 


(84) 


(86) 


(86) 


(87) 


6k  -  6  (24) 

appears  quite  reasonable  since  the  ratio  does  not  depend  upon  the  time  step  [2]. 
Factoring  out  an  “e  *T/2”  in  the  denominator,  Equation  (87)  becomes 
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aT  e  ' 

(eaT/2  -  e~ aT/2  j 

and,  since, 

aT/2 

e  -  sinh  aT/2  +  cosh  aT/2 

and 

sinh  aT/2  ■  1/2  ^e'17^2  «  e-aT/2  j 

after  a  few  manipulations  one  has 

(aT  / 2)  [ctnh  (aT/2  )  +  1] 

the  ratio  for  an  exponential  input  to  an  Eular  integrator. 
If  “a”  is  imaginary,  that  is, 

a  -  i  to  , 

Equation  (91)  becomes 

(u)T/2)  [ctn  (d)T/2)  +  i] 


(88) 


(89) 


(90) 


(91) 


(92) 


(93) 
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Since  this  “ratio”  is  complex,  one  must  determine  its  amplitude  and  phase. 
Multiplying  Equation  (93)  by  its  complex  conjugate  and  taking  the  square  root,  one 
has 


(uiT/2)  |ctn2(u)T/2)  +  1]1/2 


(94) 


which  readily  simplifies  to 


( u)T/2) 
sin(u)T/2  ) 


(95) 


for  the  amplitude  ratio  of  a  sine  wave  input  to  an  Eular  integrator. 

Dividing  the  imaginary  part  by  the  real  in  Equation  (93),  and  taking  the  arc 
tangent  yields 


(uT/2) 


(96) 


for  the  phase  error  of  a  sine  wave  input  to  a  Eular  integrator. 

One  would  proceed  in  a  similar  manner  for  the  other  integrators  of  interest;  see 
Reference  2  (6,  and  10)  for  details.  The  results  are  summarized  in  Tables  5  and  6. 

For  a  low  frequency  sine  wave  input  the  Eular  integrator  has  half  the  error  of  the 
Trapezoidal  integrator,  but  unfortunately  it  has  a  linear  phase  error.  For  6  =  1/2, 
MVC  does  not  have  the  phase  error  and  since  there  would  be  no  advantage  to 
choosing  a  “6”  other  than  one  half,  it  is  worthy  of  consideration.  Also, the  ratio  of 
MVC  error  to  TC  error  is  1:2,  and  the  errors  are  of  opposite  sign  which  serves  to 
explain  their  mixture  in  R-KC. 


The  linear  phase  error  of  EC  and/or  the  fact  that  MVC  samples  at  mid  interval 
for  6  =  1/2  may  explain  the  observation  that  in  some  simulations  the  difference 
between  the  simulation  and  exact  answer  is  a  shift  in  time. 


TABLE  5.  ACCURACY  FOR  AN  EXPONENTIAL  INPUT  TO  AN  INTEGRATOR 


T 


TABLE  6.  ACCURACY  FOR  A  SINE  WAVE  INPUT  TO  AN  INTEGRATOR 


Though  it  is  unwise  to  sample  at  the  Shannon  limit  of  two  samples  per  cycle,  it  is 
instructive  to  examine  the  ratios  and  phases  at  that  limit,  Table  8.  The  noisy 
qualities  of  the  Simpson  integrator  are  apparent.  Since  R-KC  takes  twice  as  many 
samples  as  Simpson,  it  has  no  difficulty  at  the  limit;  its  ratio  is  two  thirds  that  of  M  VC 
since  TC’s  is  zero.  Trapezoidal  integration  appears  particularly  bad  at  the  limit  since  it 
might  stabilize  an  otherwise  unstable  system.  The  last  entry  was  tuned  [2]  for  unity 
ratio  at  the  limit  but  it  has  a  ninety  degree  phase  error.  The  implications  in  determining 
gain  and  phase  margins  are  obvious. 


TABLE  8.  ACCURACY  AT  THE  SHANNON  LIMIT  ( wT  =  tt  ) 


RATIO 

PHASE 

Exact 

1 

(0  dB) 

0 

M.V.C. 

~  =  1.57... 

(2  dB) 

(1/2  -  6)tt 

E.C. 

y  -  1.57... 

(2  dB) 

TT 

~  2 

T.C. 

0 

(-  dB) 

0 

Simpson 

00 

(°°  dB) 

0 

R-KC. 

TT  _  _  _ 

•j  *  1.05... 

(0.2  dB) 

0 

E.T.C. 

(A- 6)tt 

♦  1 
-  2 

E.T.C. 

6  *  -+  f-)1^2 

2  V 

1.28 

(1  dB) 

♦  1 
'  2 

E.T.C- 

1 

(0  dB) 

♦  1 
"  2 
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6.  SOME  SAMPLE  PROBLEMS 


Now  to  the  meat  of  the  subject,  the  numerical  solution  of  differential  equations.  To 
incorporate  non-zero  initial  conditions,  one  should  proceed  from  the  differential 
equation  using 


LfX(n)(t)  ] 


n 


s 


n-1 

X(s)  -  l 

£-0 


n-£  -1 
s 


x«> 


(97) 


To  attempt  to  proceed  directly  from  the  transfer  function  is  difficult  because  to  obtain 
the  transfer  function,  the  ratio  of  output  to  input,  it  was  assuned  that  the  initial 
conditions  were  all  zero. 

A.  Single  Integration 

For  “n  =1”  Equation  (97)  becomes 


X(s)  -  s  X(s)  -  X 


(98) 


and  solving  for  “X(s)”  one  has 


X(s)  «  +  —  • 

s  s 


(99) 


laking  the  z-transform  of  Equation  (99), 


X(z)  -  z(^)  +  XoZ^). 


(100) 


There  is  no  difficulty  in  taking  the  z-transform  of  the  initial  conditions  since  they 
are  constants  in  the  frequency  domain;  they  are  multiplied  by  Dirac  deltas  in  the 
time  domain  [2]. 
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The  input,  “XU)",  is  another  matter  and  requires  the  application  of  the 
convolution  approximations,  Table  2. 

Now  to  develop  the  recurrence  for  MVC,  Table  9, 

(l-z)X(z)  =  T  z  X(z,l/2)  +  Xo  (101) 


TABLE  9.  CONVOLUTION  APPROXIMATIONS  FOR  SINGLE  INTEGRATION 


z  («‘)\ 

M.V.C. 

t  x(z,  1/2)  r  1  1 

L1  - 2  J 

E.C. 

<lIz)lX<2>  '  V  °rTX(2>  [l-,’1] 

E.T.C. 

Tin  +  (l-n)zl  x(z)  T  Xo 

1 -z  1-z 

T.C. 

R-KC. 

6(l-z)  [z*(z)  +  4  (XCz>  '  1/2)"  x<-T/2.)j+  (X(z)  -  Xo)] 

and  substituting  the  definition  of  the  (modified)  z-transform,  one  has 

CD  00 

l  X(nT)  zn-  X(nT)  zn+1  =  t£  &(nT+y )  zn+1  +  X  •  (102) 

n=0  n=0  1  ° 


Now,  the  important  step,  equating  coefficients  of  like  powers  of  “z",  one  has 


X0  *  Xo 


(103) 


X  *X  ,  +TX  ...  ,  n>0 

n  n-1  n-1/2 


(104) 


It  is  important  to  note  that  the  general  recurrence,  Equation  (104),  does  not  apply  at 
“n  =  0". 

For  EC  there  are  two  possibilities  because  of  the  lack  of  symmetry. 


(1-z)  X(z)  =  T  (X(z)-*  ]  +  X 

o  o 


(105) 


becomes 


l  X(nT)zn  -X(nT)  zn+1 =  T £  X (nT)  zn  + 


n-0 


n=l 


(106) 


and  equating  coefficients  of  like  powers  of  “z". 


x  =  x 

o  o 


(107) 


*„  ■  Vi +  T  *„•  ">0  • 


(108) 


Rectangular  Integration! 


For 


(1-z)  X(z)  -  T  z  X(z)  +  Xq 


(109) 


36 


(110) 


}'  X(nT )  zn  -  X(nl)  z"+1  =  T  )  X(nT)  z°+1  +  X0 
n=0  n=0 


and . finally. 


x 


o 


dll) 


X 

n 


X  ,  +  T  X  ,  n>0 

n-1  n-1. 


(112) 


\Eular  Integration. 

\ 

Proceeding  in  a  similar  manner  would  produce  Table  10.  It  might  appear  at  first 
that  MVC  and  R-KC  are  not  symmetrical  since,  if  the  rolls  of  “g(s)”  and  “f(s)”  are 
interchanged,  one  would  have 


[X(z ,  -  1/2)  -  X(T/2)  ] 

I  -  z 


(113) 


and 


6  J_2)  [  (*(z)  -  XG)  +  4z  X(z,l/2)  +  z  X(z)  ] 


(114) 


respectively,  but  these  approximations  lead  to  the  same  recurrences  found  in 
Table  '  9. 

Of  course,  these  recurrences.  Table  10,  are  those  used  to  develop  the  convolution 
approximations  in  an  earlier  section,  and  are  not  unexpected.  At  least  they 
demonstrate  consistency. 
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TABLE  10.  RECURRENCES  FOR  SINGLE  INTEGRATION 


Xo  =  xo 

M.V.C. 

X  =  X  ,  +  T  X  1  /0 
n  n-1  n-1/ 2 

E.C. 

X  *  x„  ,  +  T  X„ 

n  n-1  n 

E.T.C. 

x  =  x  ,+T[nx  +  (i-n)  x  .  ] 

n  n-1  n  n-1 

T.C. 

X»  ■  Xn-1  +  R  +  Vl> 

T.T.C 

Xn  =  Xn-1  +  Ifc-W  V  <1+2r»  Xn-1+  2Vl+  Xn-2)] 

R-KC. 

X„  -  Vl  +  !  +  «n-l/2  +  Vl> 

B.  Double  Integration 

Double  integration  is  interesting  because  it  illustrates  some  of  the  problems 
associated  with  initial  conditions  and  start  up.  Equation  (97)  for  “n  =  2”  is 

X(s)  =  s2X(s)  -  s  Xq  -  Xq  (115) 

and  solving  for  "X(s)”,  one  has 

X(s)  -  +  if  +  ^  .  (116) 

i  s 
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Taking  the  z-transform  of  Equation  (116), 


si"'  -  f  fi +  Ml-H  z(»)- 


(117) 


Before  proceeding,  consider  the  case  where  )t(t)  =  0,  that  is , 


«-•>  -  7(.!')  *  \>  7(» 


(118) 


The  exact  z-transform  is 


TzX  X 

XU>  -  +  (T-z)  * 

(i-zr  1  ’ 


(119) 


or,  rearranging. 


(I  -  z)~  X(z)  -  T  z  X  +  (l  -  z)Xo 


(120) 


Substituting  the  definition  of  the  z-transform  and  equating  coefficients  of  like 
powers  of  "z”, 


*o  *  V 


(121) 


X.  -  X„  +  T  X 

1  O  D 


(122) 


X  -  2  X  .  -X  ...  n>l. 
n  n-I  n-2 


(123) 


It  is  important  to  note  that  the  general  recurrence,  Equation  (123),  may  be  first 
applied  at  “n  =  2”  and  DOES  NOT  apply  at  “n  =  1”.  As  the  order  of  the  system 
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1 


increases  the  number  of  steps  before  the  general  recurrence  may  be  applied 
increases!  For  further  discussion  on  this  point  see  Reference  2,  in  particular, 
Appendix  D. 

With  a  little  mathematical  induction  one  may  readily  demonstrate  that 


*n 


X  +  n  T  X  iiM) 
o  o. 


(124) 


Mathematical  induction  serves  to  further  illustrate  the  point  made  above.  In 
addition  to  showing  a  relationship  for  *‘n"  is  valid  for  “n  +  1”,  one  must  also  find 
some“n”  for  which  it  is  valid.  For  equation  (124)  the  smallest  such  “n”is“n=  1”  but 
for  Equation  (123)  it  is  “n  =  2”. 

Table  11  contains  the  required  approximations  convolution  and  Table  12  the 
recurrences.  Note  how  the  input  enters  the  startup  step,  “n=l”  in  particular  for  TC. 

TABLE  11.  CONVOLUTION  APPROXIMATIONS  FOR  DOUBLE  INTEGRATION 
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7 

•te*. 


TABLE  12.  RECURRENCES  FOR  DOUBLE  INTEGRATION 


C.  Single  Pole  Filter 

Now  for  a  bone  with  a  little  more  meat  on  it,  the  single  pole  filter  whose 
differential  equation  is 

X(t)  -  -a  X(t)  +  g(t).  (125) 

Applying  Equation  (97),  one  has 

s  X(s)  -  XQ  *=  -a  X(s)  +  g(s),  (126) 
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and  solving  for  “X(s 


X(s) 


g(s)  +  X 

_ O 

S  +  ti 


(127) 


For  X0  =  0,  one  would  have  the  transfer  function 


X(s)  _  1_ 

g(s)  s  +  a 


(128) 


taking  the  z-tranaform  of  Equation  (127), 

«z)  -  z(!&>) + 

Proceeding  as  before,  one  has  Tables  13  and  14.  Like  the  single  integrator  there  is 
no  startup  step  but  note  how  past  values  are  “decayed”  based  on  how  "old”  they 
are. 

A  lead-lag,  in  Equation  (130),  is  sometimes  required. 


(130) 


To  incorporate  initial  conditions  convert  Equation  (130)  back  to  the  time  domain. 


X(t)  +  a  X(t)  »  g(t)  +  b  g(t). 


(131) 


and  then  apply  Equation  (97), 


TABLE  13.  CONVOLUTION  APPROXIMATIONS  FOR  A  SINGLE  POLE  FILTER 


M.V.C. 

T  -*T/2 

— - —  <g(z.  ”1/2)  -  g(-T/2) ) 

1-z  e 

E.C. 

T  z  e  aTg ( z ) 

(1  -z  e_aT) 

T.C. 

|(l  +  ze_aT) 

_aT  ^  ~  80] 

(1  -  z  e  ) 

_aT 

R-KC. 

T  vr  1(1+2  e  aT)  8(z>  ~  80  +  4e  (g(z,  -1/2)  -g(-l/2T))] 

,6(1  -z  e) 

TABLE  14.  RECURRENCES  FOR  A  SINGLE  POLE  FILTER 
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and  solving  Equation  (132)  for  “X(s)’ 


(s  +  b)  g(s)  -  g  +  X 

X(s)  -  - - - 2 - £ 

s  +  a 


(133) 


taking  the  z-transform  of  Equation  (133), 


x(z>  -  z[$^*<s>] +  <V«o>  ?(—) 


(134) 


Noting  that 


s  +  b  j  (a  -  b)  , 

s  +  b  s  +  a 


(135) 


Equation  (134)  may  be  written 


X(z)  •  g(z)-Z[^|g(s)]  +  <Xo  -  «„)  Z^)  - 


(136) 


The  approximation  required  in  Equation  (136)  is  the  same  as  the  single  pole  filter 
with  the  coefficient  “(a-b)”.  Table  15  contains  the  recurrences. 

If  the  lead-lag  is  -of  the  form, 


1  +  s/b  m  a  /s  +  b\ 
1  +  s/a  b  \s  +  a/ 


(137) 


then  the  recurrences  in  Table  16  would  apply. 


If  one  had  the  form 


8 

1  +  s/a 


(138) 
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in  Table  16  let 


(a  -  b)  ■+  a  •  (139) 


(140) 


Unfortunately  Table  17  does  not  lend  itself  to  differentiation  since  the  requirement 
that  “aT«l”  results  in  slow  response.  (See  Appendix  C  for  some  comments  on 
differentiation.) 


TABLE  17.  RECURRENCE  FOR  /  -S-2- 

\*  +  « 


M.V.C. 

Xn  =  a(V  aT  e"aT/2gn-l/2)+e"al(Xn-l  "  a  gn-l> 

E.C. 

x„-»(l-aT)  V  ="T<Vr**,-,> 

T.C. 

Xn  “  a  (1-  aT/2  )  8n+e‘aT(Xn_l-a(1  +  aT'/2) 

R-KC. 

Xn  -  .(1  +.1/6)  g„-iaij£:aT/2  gn_1/2 

+e  ~aT  <Vi  -  a  (1  +  aT/6)  gn-1 ) 
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D.  A  Forced  Damped  Oscillator 

This  problem  was  discussed  in  some  detail  in  Reference  2,  and  only  some  results 
will  be  presented  here.  The  transfer  function 


s2  +  2  C  uos  +  «* 


,  1  >  c  >  0 


is  represented  as 


2  2 
(s  +  a)  +0J 


where 


(141) 


(142) 


a  -  C<jJ 


(143) 


and 


0)  = 


ci-c  V/2« 


(144) 


the  startup  step(s)  will  be  found  in  Table  18  and  the  general  recurrence^)  in  Table 
19. 

Figure  2  serves  to  illustrate  the  difference  between  the  innerconnected  integrator 
approach,  Figure  2(a),  and  that  suggested  herein,  Figure  2(b).  The  velocity,  “X(s)”, 
would  be  found  using  the  single  pole  filter  recurrences,  Table  14.  Unless  “£(s)”  and 
“X(s)”  are  needed  for  some  other  computation  or  output,  they  need  not  be  computed. 
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TABLE  18.  STARTUP  STEP(S)  FOR  OSCILLATOR 


X  -  X 
o  o 


M.V.C. 


E.C. 


T.C . 


R-KC 


X.  -  e_aT(cos  U)T  +  ^sin  toT)  X  +  e  _aT  .sJLn,M  X 
1  03  o  0)  C 


.  _  -  aT/2  sin  oiT/2 
+  l  e  - q 

0)  i/2 


+  0 


T  -  aT 

+  2  6 


sin  o)T 

0) 


+  ^  e'aT/2  sin  ojT/2  ...  -aT/2  , 

3  - -  [2  g  ,  +  e  cos  u>T/2  q  ] 

u)  vl/2  0 


TABLE  19.  GENERAL  RECURRENCES  FOR  OSCILLATOR,  n>1 


M.V.C. 


E.C. 


T.C. 


R-KC. 


_  -aT  _  v  -2aT  _ 

X  M  2c  cos  uT  X  . —  c  X  m  .... 
n  n-i  n— i 


T*  «  [*n-l/2+  *  *n-3/2j 


.aT 

+  T  e  sin  oiT  g  , 
-  n— 1 


.  _  -aT  sin  u)T  . 

+  T*  — 3 -  ®n-l 


2T  -«T/2  sln(  hff/2  )  /  -«T/2  -«T  \ 

+  “6  -  (*„-l/2+e  '  cos  utf/2  In-i^  ‘T  V3/2| 
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Figure  2a.  Interconnected  tingle  Integrator! 


Figure  2.  A  forced  damped  oscillator. 


7.  RESULTS  AND  CONCLUSIONS 


The  Heaviside  unit  step  and  the  Dirac  delta  were  shown  to  be  useful  in  outlining 
a  proof  of  the  Ragazzini-Zadeh  identity.  The  mean  value  theorem  of  the  integral 
calculus  and  various  numerical  integrators  were  used  to  develop  approximations 
of  use  in  the  digital  simulation  of  continuous  systems.  The  accuracy  of  these 
approximations  was  then  studied  for  various  inputs  to  a  single  integrator.  The 
inputs  considered  were  a  constant,  ramp,  parabolic,  exponential  and  sine  wave. 
Recurrences  were  then  developed  for  single  and  double  integration  to  illustrate 
incorporation  of  initial  conditions.  Recurrences  for  a  lead-lag  transfer  function 
and  a  forced  damped  oscillator  were  also  presented. 

The  sample  problems  were  simple  to  illustrate  the  approach.  The  approach  has 
been  applied  to  the  real  time  digital  portion  of  the  hybrid  simulation  of  a  high  spin 
rate  air  defense  missile  [11],  The  simulation  is  presently  in  operation. 

Another  application  would  be  subroutines  of  often  countered  transfer  functions 
for  use  in  simulation  models  such  as  Reference  12. 

Further  effort  on  the  recommendations  [2]  is  required.  It  is  hoped  that  in  trying  to 
determine  the  facts  that  additional  fallacies  have  not  been  introduced.  Comments 
are  welcomed. 


APPENDIX  A. 


TRANSFORM  TABLE  FOR  SELECTED 
FUNCTIONS*  AND  DISTRIBUTIONS 


•Mealy.  M.,  Tables  of  Laplace.  Fourier  and  Z-tramforms.  London:  W.  ind  R.  Chamber*  Ltd.,  1967.  Cadzow,  J.A.,  Discrete- 
Time  Systems.  Englewood  Cliff*,  New  Jeraey:  Prentice  Halt.  I97J. 


f(t)  f(8)  f(Z> 


r 


APPENDIX  B. 


DIRAC  DELTA  DISTRIBUTION 


55 


APPENDIX  B.  DIRAC  DELTA  DISTRIBUTION* 


CO 

r 

1  f(t)  6(t  -  nT)dt  -  f (nT) 

J 

(B-l) 

•CO 

where  f(e)  is  a  well-defined  function  at  t  »  nT. 

B(t)  -  6(-t) 

(B-2) 

6<at) "  ul 8<t) 

(B-3) 

b[8^c>]  “2lg'(nT)|  “  nl)'  g' (nT)  *  0 

n 

(B-A) 

t  8(t)  -  0 

(B-5) 

f(t)6(t  -  nT)  -  f(nT)6(t  -  nT) 

(B-6) 

Jb( t  -  x)6(t  -  nT)dt  -  6 (t  -  nT) 

(B-7) 

f  <m>  m  <m) 

I  6(t)  f(t)dt-  (-1)  f(o) 

J 

(B-8) 

-OO 

F  <l> 

/  6(o)  e”8  dt  •  s 

(B-9) 

L 

•Dirac,  P.  A.  M.,  The  Principles  of  Quantum  Mechanics,  London:  Oxford  University  Press, 
1968,  Messiah,  A.,  Quantum  Mechanics,  New  York:  John  Wiley  and  Sons,  1964. 


Differentiation  is  to  be  avoided  on  analog  and  digital  computers 
but  there  are  occasions  . . . 

Equation  (97)  for  “n  =  1”  is 

X(s)  =  s  X(s)  -  X(0) 

and  taking  the  z-transforra 

X(z)  =  Z(s  X(s) )  -  X(0) 

Since 

-St 

z  =  e 

it  follows  that 


and  Equation  (C-2)  becomes 

X(z)  -  Z  Cn(l/z)  X(s)^  -  X(0) 

Applying  the  Ragazzini-Zadeh  identity.  Equation  (C-5)  becomes 

X(z)  -  |  inLlzj  X(z)  -  X(o) 


whenever  possible, 


(C-l) 


(C-2) 


(C-3) 


(C-4) 


(C-5) 


(C-6) 
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One  of  the  possible  series  for  the  logarithm  is 


Cn 


(1  -  z) 


n=l 


(C-7) 


and  substituting  into  Equation  (C-6)  and  equating  coefficients  of  like  powers  of  z : 


X  =  o  , 
o 

X.  -  X 

X  =  -i _ £ 

*1  T 


x7  -  X.  (X7  -  X.)  -  (X.  -  X  ) 
y  _  1  1  .  I  1 _ 1  O  » 

2  “  T  2T 


X3  -  X2  (X3  -  X2)  -  (X2  -  Xx) 
X3  =  T  +  2T 


+  [(X3  -  x2)  -  (X2  -  X3)]  -  [<X2  -  xx)  -  (X3  -  xo)] 


3T 


(C-8) 

(C-9) 

(C-10) 


(C-ll) 


This  series  was  chosen  because  each  term  is  a  finite  difference,  and  there  is  no  attempt 
to  use  any  term  unless  the  whole  term  can  be  used.  That  is,  for  a  given  “Xn," 


tn 


(C-12) 


Since  all  past  values  of  “Xn”  are  used,  a  general  recurrence  never  occurs! 
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One  might  truncate  the  series  at  one,  two  or  three  terms,  i.e.. 


n  >  0 


(C-13) 


X  = 
n 


3  X  -  4  X„  .  +  X 
n _ n-1 _ n- J 

2T 


,  n  >  1 


(C-14) 


11  X  -  18  X  ,  +  9  X^  „  -  2  X  , 
n _ n-1 _ n-2 _ n-J 

6T 


n  >  2 


(015) 


respectively. 

Of  course,  integrators  may  be  developed  in  a  similar  fashion  (Table  C-l). 


TABLE  C-1  INTEGRATORS  VIA  LOG  APPROXIMATION 


£n(l/z) 


—  =  T  fn(Vz)"1 
s 


CLASSICAL 

NAME 


1  -  z 


1  -  z 


irfi) 


T  z 
1  -  z 


1  -  z 


I  (Hi) 


3T  (1  +  z)' 
8(1  -  z3) 


Eular 


Rectangular 


Trapezoidal 


Simpson's  3/8 


Since  some  checks  are  warranted,  consider  the  differentiation  of  a  cosine  wave. 
Equation  (C-13)  would  become 


cos  nu)T  -  cos  (n  -l)uiT 

T  . . 


(C-16) 


which  may  be  reduced  to 


-u>  ( S' 1  /  2 2  — )s  1  n  < moT  -  U)T/2  ) 


(C-17) 


Of  course  the  answer  is  “ho  sin(n<oT)’\ 


The  mean  value  theorem  of  the  differential  calculus  states 


>  -  *n  *n-l  ,  0  <_  n  <_  1. 

n-n  nT  -  (n-l)T 


(C-18) 


If  in  Equation  (C-18) 


n  =  1/2 


(C-19) 


then  in  Equation  (C-17) 


s'n  mT/2 
u>T/2 


(C-20) 


may  be  interpreted  as  an  amplitude  error. 
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Since  there  is  tunable  integration  [5],  why  not  tunable  differentiation? 


Y(1  -z  ) 
Tz 


(1  -  Y)U  -  z)  _ 
T 


Y  +  (1  -  2y)z  -  (1  -  y) z 
Tz 


Invoking  the  mean  value  theorem,  for  a  cosine  wave 


-to  sin  (nu)T  -  n^T)  =  y  cos  (n  +  l)u)T  +  (1  -  2y)  cos  nu)T 


-  (1  -  Y)  cos  (n  -  1)ojT 


sin  noiT 


and  finally 


tan  nwT  =  (2y  -  1)  tan  wT/2 


(C-21) 


(C-22) 


(C-23) 


(C-24) 


Therefore 


or 


> 


1  /  2 1 1  + 


tan 


tan 


tyr _ 

u)T/2 


(C-26) 


In  small  angles,  “uT  -  n 


n  '  'i  -  1/2 


(C-27) 


or 


>  '  1/2  +  n  (C-28) 

For  any  “wT,” 


TABLE  C-2  TUNINGS  INDEPENDENT  OP  ‘  wT” 


n 

Y 

-  1/2 

o  ! 

0 

1/2 

1/2 

1 

When  the  other  differentiators  of  interest  are  checked  in  a  manner  similar  to 
Equation  (C-I3),  Tables  (C-3)  and  (C-4)  would  rest  It.  The  only  difference  between 
the  tables  is  one  of  point  of  view.  When  small  angle  approximations  are  warranted, 
Tables  (C-5)  and  (C-6)  are  appropriate. 

One  may  conclude  that  for  the  first  two  differentiators,  “y  =  0”  and  “y  =  I,"  a 
phase  shifted  interpretation  is  more  appropriate,  while  for  the  last  two,  *‘n  =  2"  and 
“n  =  3,”  it  is  not.  At  the  Shannon  limit,  Table  (C-7),  the  first  two  are  36%  low  while 
the  last  two  have  suffered  a  serious  “phase”  error.  Differentiation  is  to  be  avoided. 


DIFFERENTIATION  OF  A  COSINE  WAVE,  “AMPLITUDE/  PHASE 


sin  (m  'T  +  3  wT/ 2  ) 


TABLE  C-4  DIFFERENTIATION  OF  A  COSINE  WAVE 


cos  nwT 


TABLE  C-5  DIFFERENTIATION  OF  A  COSINE  WAVE,  uT  - 


Z/l"t  +  I"u)  uiS| 


TABLE  C4  DIFFERENTIATION  OF  A  COSINE  WAVE,  u>T  '  * 
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S  t'—ICJ  3  - - - 


TABLE  C-7  DIFFERENTIATION  OF  A  COSINE  WAVE,  u7  -  jt 


f 

-a)  sin  nu)T 

1  -  z 

T 

(Y  =  0) 

-0)(2/n)  sin(na)T  -tt/2) 

1  -  z 

Tz 

(Y  =  1) 

-uj(2/n)  sin(nuiT  +  71/2 

1  -  z2 
2Tz 

(Y  *  1/2) 

0 

z~1/2  - 

1/2 

z 

-uj(2/77  )  sin  nu>T 

T 

I  f  XL 

T  L 
n=l 

7 Z>  <n  "  2> 

-u)  (-4/ti  )  cos  nu)T 

I  y  XL 

n=l 

7^  <--» 

-o)(  20/377  )  cos  noiT 

APPENDIX  D.  R-K(N)  CONVOLUTION 


Note  the  rather  poor  behavior  at  wT  =  2tt  /  3  ;  BOOM! 


One  of  the  forms  used  in  fourth  order  Runge-Kutta  integration 


and  the  ratio  for  a  sine  wave  input  is 


1  (a£\ 

4  \  2  / 


cos  a)T/2  +3  cos  u)T/  6 


sin  uiT/2 


At  the  Shannon  limit.  “u)T  -  tt  the  ratio  is  1.0202621 . . . 


Equation  (D-3)  may  be  written 


X(z) 


(, 

Ft  /i  +  z\" 

,  3 

~Tz1/3' 

.  3 

rTz2/3ii 

I4 

[2  \1  ~  *) 

+  8 

1-2 

+  8 

Li  -  Jl 

X(z)  + 


X(0) 
1  -  z 


(D-5) 


As  with  R-K(3),  Equation  (D-5)  is  a  linear  combination  of  trapezoidal  integration 
and  mean  value  integration  but  with  7  =  1/3,  2/3  in  this  case.  Therefore 


’jg(s)  f (s) J  ~  |  |g(z) 


(f(z)  -  f  (0) )  +  3  g(zM)  (f  (z  ,-1/3)  -  f  (-T/3)) 


+  3  g(z,  2/3)  (f  (z  ,  -2/3)  -  f(-2T/3)) 


+  f(z)  (g(z)  -  g(0))  ,  (D-6) 


R-K(4)  Convolution. 


Higher  order  convolutions,  R-K(N)C,  could  be  developed  in  this  fashion  using 
trapezoidal  integration/convolution  and  mean  va’ue  integration/convolution. 
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Eglin  Air  Force  Bate.  Florida  32542 

DRSMI-LP.  Mr  Voigt  1 

DRDMI-T.  Dr  Kobler  1 

DRDMI-TE  1 

DRDMI-TG,  D.  Clliax  1 

DROMI-TD  1 

DRDMI-TDD  1 

DRDMI-TDK  1 

DROMI-TDS  1 

DRDMI-TDR  1 

DRDMI-TDF.  V  Grimes  1 

DRDMI-TDD.  R.  Dickson  100 

DRDMI-Tl  1 

DRDMI-TK  1 

DRDMI-TB  5 

DRDMI-TBD  3 

DRDMI-Tl  (Reference  Copy)  1 

DRDMI-Tl  (Record  Copy)  1 
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